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ABSTRACT 

A  general  description  is  provided  of  current  model  development  associ¬ 
ated  with  Task  DST  98/091  Modelling  Airbreathing  Combustion.  This  note 
describes  current  and  future  work  towards  the  completion  of  a  fully  imple¬ 
mented  software  suite  for  turbulent  combustor  applications.  Particular  at¬ 
tention  is  paid  to  the  issues  surrounding  applying  the  Conditional  Moment 
Closure  (CMC)  family  of  models,  which  were  developed  in  simple  flows,  to 
flow  cases  of  practical  interest.  These  issues  include  the  integration  of  radi¬ 
ation  heat  transfer,  conductive  heat  transfer,  complex  recirculating  flow,  and 
complex  chemistry  within  the  CMC  model  context. 
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Issues  Surrounding  the  Development  of  Models  for  use  in 
Predicting  Gas  Turbine  Combustion 


EXECUTIVE  SUMMARY 

The  prediction  of  combustion  within  gas  turbine  engine  combustors  (GTCs)  is  a  criti¬ 
cal  requirement  for  understanding  how  to  successfully  operate  gas  turbine  engines  and  the 
platforms  in  which  they  are  installed.  The  combustion  processes  which  occur  in  GTCs  gov¬ 
ern  important  operational  issues  such  as  engine  component  lifetimes,  engine  performance, 
stability  of  engine  operation,  and  detectable  emissions.  Thus  the  ability  to  comprehen¬ 
sively  predict  all  aspects  of  gas  turbine  combustion  can  usefully  aid  in  the  diagnoses  of 
poor  engine  performance  and  unusual  component  degradation,  as  well  as  provide  inputs 
to  the  estimation  of  platform  infra-red  and  visible  signatures. 

By  developing  a  modelling  philosophy  based  upon  fundamental  physical  processes, 
rather  than  empirical  observation  of  particular  engines,  an  inherent  flexibility  is  lent  to 
the  final  model.  This  allows  it  be  to  applied  across  all  existing  engine  types,  and  indeed,  to 
predictions  of  the  characteristics  of  foreign-operated  and  future  engines.  Task  DST  98/091 
Modelling  Airbreathing  Combustion  was  established  to  provide  DSTO  and  the  ADO  with 
a  fundamentally  based  GTC  predictive  capability,  implemented  as  a  suite  of  numerical 
software.  After  one  year  of  development,  it  is  appropriate  to  pause  and  delineate  the  way 
forward  in  the  research  program.  The  purpose  of  this  technical  note  is  to  provide  just  such 
a  technical  delineation. 

The  key  physical  processes  which  contribute  to  modelling  GTCs  are  outlined  in  this 
note.  These  processes  are  described  to  result  from  interactions  between  four  key  elements 
of  the  overall  physical  system.  These  key  elements,  turbulence,  combustion,  heat  transfer, 
and  multi-phase  flow,  are  each  relatively  well  understood  in  isolation.  While  a  myriad  of 
mature  modelling  methodologies  are  available  for  each  element  in  isolation,  the  integration 
of  different  models  from  different  development  backgrounds  takes  some  care.  This  note 
describes  in  detail,  some  of  the  issues  surrounding  the  implementation  of  a  number  of 
disparate  models  in  a  coherent  model  framework  capable  of  being  applied  to  gas  turbine 
combustion  prediction. 
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1  Introduction 

The  prediction  of  combustion  within  gas  turbine  engine  combustors  (GTCs)  is  a  criti¬ 
cal  requirement  for  understanding  how  to  successfully  operate  gas  turbine  engines  and  the 
platforms  in  which  they  are  installed.  The  combustion  processes  which  occur  in  GTCs  gov¬ 
ern  important  operational  issues  such  as  engine  component  lifetimes,  engine  performance, 
stability  of  engine  operation,  and  detectable  emissions.  Thus  the  ability  to  comprehen¬ 
sively  predict  all  aspects  of  gas  turbine  combustion  can  usefully  aid  in  the  diagnoses  of 
poor  engine  performance  and  unusual  component  degradation,  as  well  as  provide  inputs 
to  the  estimation  of  platform  infra-red  and  visible  signatures. 

By  developing  a  modelling  philosophy  based  upon  fundamental  physical  processes, 
rather  than  empirical  observation  of  particular  engines,  an  inherent  flexibility  is  lent  to 
the  final  model.  This  allows  it  be  to  applied  across  all  existing  engine  types,  and  indeed, 
to  predictions  of  the  characteristics  of  foreign-operated  and  future  engines. 

Currently,  DSTO  has  a  Task  in  place  to  develop  a  software  suite  to  numerically  predict 
combustion  characteristics  within  GTCs  and  other  steady  constant-pressure  nonpremixed 
combustion  devices.  Task  DST  98/091  Modelling  Airbreathing  Combustion  is  at  the  end 
of  its  first  year  of  execution.  This  note  is  designed  to  broadly  outline  the  theoretical 
and  technical  issues  which  surround  the  continuing  effort  to  develop  and  implement  GTC 
model  software. 

The  GTCs  currently  employed  in  modern  aeropropulsion  applications  employ  a  non¬ 
premixed  combustion  configuration.  In  this  configuration,  liquid  fuel  is  sprayed  into  a 
stabilised  recirculating  air  flow.  The  liquid  fuel  evaporates,  mixes  with  surrounding  hot 
gases,  and  burns  concurrently  within  the  fuel-rich  recirculation  zone.  This  type  of  config¬ 
uration  is  employed  in  preference  to  premixed  fuel-air  delivery  due  to  its  inherent  stability 
of  operation  over  a  wide  range  of  power  settings.  The  combination  of  a  low-speed  recir¬ 
culation  zone  and  a  localised  abundance  of  fuel  in  the  zone  ensure  the  presence  of  a  flame 
under  these  various  operating  conditions.  Further  air  is  added  downstream  of  the  recir¬ 
culation  zone  to  complete  the  combustion  of  fuel  present  in  the  product  gas  stream.  Still 
further  air  is  progressively  added  to  cool  the  stream  to  a  point  where  the  gas  temperature 
is  low  enough  to  avoid  damage  to  downstream  components  in  the  turbine  assembly. 

The  overall  process  of  heat  injection  into  the  stream  passing  through  the  combustor 
involves  many  interacting  physical  processes.  Intense  turbulent  mixing  is  required  within 
the  combustor  to  bring  fuel-rich  and  fuel-lean  pockets  of  gas  together  at  a  rapid  rate.  This 
allows  combustion  reactions  to  proceed  and  the  rapid  release  of  heat  in  a  confined  volume 
to  take  place.  The  overall  chemical  reaction  which  is  taking  place  within  the  combustor 
is  only  a  global  representation  of  the  multitude  of  elementary  reactions  going  on  between 
many  thousands  of  intermediate  chemical  species.  These  reactions  proceed  at  rates  which 
are  non-linearly  dependent  upon  temperature  and  species  concentrations,  and  are  variously 
accelerated  and  impeded  by  the  turbulent  mixing  processes  which  continously  change  the 
distributions  of  these  quantities.  To  further  complicate  the  picture,  the  hot  combusting 
mixture  is  continuously  exchanging  heat  with  the  surrounding  combustor  surfaces  through 
radiative  absorption.  The  formation  of  soot  particles  from  fuel  rich  gases  significantly 
enhances  this  process,  and  at  the  same  time  adds  the  complication  of  a  reacting  solid- 
phase  to  the  problem  to  be  modelled.  Transfer  of  heat  to  liquid  fuel  droplets  and  their 
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evaporation  adds  another  aspect  to  the  overall  process. 


Though  there  are  many  different  physical  processes  inherent  in  GTC  combustion,  this 
note  is  limited  in  scope  to  separate  discussions  of  a  selected  number  of  the  modelling 
methods  which  will  be  implemented  in  the  near  future. 


2  Modelling  Methods 


Many  different  submodels  are  required  to  be  integrated  into  a  whole  in  order  to  pro¬ 
vide  a  model  framework  which  is  applicable  to  all  aspects  of  nonpremixed  airbreathing 
combustion.  In  Figure  1,  the  key  physical  elements  of  this  framework,  and  the  causal 
relationships  between  these  elements  has  been  conceptualised  according  to  the  diagram 
shown. 


PHYSICAL  ELEMENTS  OF  GAS  TURBINE  COMBUSTION 


TURBULENCE 
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Figure  1:  Schematic  diagram  indicating  the  interactions  between  key  physical  processes 
which  require  modelling  in  the  context  of  a  gas  turbine  combustor. 
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The  model  problem  has  been  broken  up  into  four  key  elements,  namely  turbulent  flow 
and  mixing  in  a  complex  geometry,  gaseous  nonpremixed  combustion,  heat  transfer  to 
boundaries,  and  multi-phase  species  dynamics.  In  Fig.  1,  the  processes  governed  by  these 
key  elements  interact  with  one  another  to  varying  degrees.  For  example  (see  horizontal 
arrow  at  top  of  the  figure),  gaseous  nonpremixed  combustion  is  influenced  by  turbulent  flow 
and  mixing.  The  range  of  scales  of  turbulent  motion  cause  pockets  of  reactive  species  to 
be  mixed  together  at  rates  which  can  limit  the  overall  rate  of  reaction  and/or  modify  flame 
structures.  Conversely,  the  heat  release  associated  with  combustion  causes  large  spatial 
variations  in  fluid  densities  which  then  affect  the  momentum  transfer  characteristics  within 
the  flow  field. 

Another  interaction  is  that  where  the  heat  release  from  gaseous  combustion  causes 
local  temperature  rises  which  lead  to  increased  radiant  and  convective  heat  transfer.  The 
presence  of  this  heat  transfer  alters  local  temperatures  which  impacts  upon  the  rates  of 
chemical  reactions  in  these  regions. 

The  presence  of  condensed-phase  particles,  such  as  soot  and  fuel  droplets,  can  cause  gas 
phase  reactions  in  their  vicinity  to  be  modified.  In  the  case  of  soot,  the  solid  phase  particles 
are  a  source  and  sink  of  chemical  species  and  enthalpy.  In  the  case  of  fuel  droplets,  the 
liquid  fuel  is  a  sink  for  sensible  enthalpy  and  a  source  of  combustible  fuel  which  evaporates 
off  the  droplet  surface  as  a  result  of  elevated  temperatures  in  the  adjacent  gas. 

The  key  to  building  the  required  model  framework  is  in  the  selection  of  compatible 
submodels  which  have  approximately  the  same  degree  of  intrinsic  accuracy  in  their  selected 
application.  The  submodel  with  the  poorest  accuracy  will  adversely  affect  the  overall 
predictive  capability,  and  so  it  is  important  to  avoid  wasting  effort  in  strengthening  strong 
links  in  the  chain  while  overlooking  weaker  links. 

The  traditionally  weakest  link  in  the  model  chain  has  been  the  modelling  of  the  gaseous 
nonpremixed  combustion  element  in  the  presence  of  turbulence.  It  is  logical  therefore  to 
start  the  assembly  of  the  model  framework  with  the  selection  of  an  accurate  turbulent 
combustion  model.  The  bulk  of  this  note  is  concerned  with  the  practical  details  of  this 
model. 


2.1  Combustion  modelling  :  CMC 

When  modelling  turbulent  nonpremixed  combustion,  difficulty  is  encountered  in  cor¬ 
rectly  determining  mean  chemical  reaction  rates.  Unlike  in  chemically  inert  turbulent 
flow,  simple  Reynolds  (or  Favre)  averaging  of  the  flow  field  does  not  provide  sufficient 
information  to  effect  an  adequate  closure  of  the  resultant  averaged  scalar  equations.  The 
highly  non-linear  dependence  of  chemical  source  terms  on  local  temperature  and  reactant 
species  concentrations,  which  fluctuate  rapidly  in  turbulent  flow,  defeats  any  attempt  at 
a  linear  first  order  closure  in  terms  of  the  averaged  local  temperature  and  reactant  con¬ 
centrations.  Many  models  for  turbulent  nonpremixed  combustion  have  been  devised  over 
the  past  fifty  years.  One  model  is  known  as  the  Conditional  Moment  Closure  (CMC) 
method  [1,  2]  and  has  been  applied  extensively  to  predicting  species  concentrations  and 
temperatures  in  gas-fuelled  turbulent  jet  flames  [3]-[8].  Instead  of  employing  conventional 
unconditional  averaging  of  the  instantaneous  reactive  species  equations,  the  CMC  method 
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involves  solving  conditionally  averaged  equations.  The  following  provides  a  description  of 
this  method,  its  rationale,  and  applications. 

The  local  instantaneous  equation  for  the  conservation  of  the  mass  fraction  ( Ya )  of  a 
chemical  species  (a)  can  be  written  as, 

I  +  £  ton) -£(*>.  57)-**  w 

where  p  is  the  fluid  density,  U{  is  the  i-th  component  of  fluid  velocity,  Va  is  the  molecular 
diffusivity  (Fickian  diffusion  approximation)  of  the  species  a,  and  wQ  is  the  net  chemical 
production  rate  of  a. 

Unconditional  averaging  of  equation  1  results  in, 

jt  «^»  +  £  <tov.»  -  ^  =  <*■>  (2) 

with  the  resultant  closure  problem  of  the  unconditional  mean  chemical  source  term, 

(pua  (Yi, . . . ,  Yn,T))  #  (p) ua  «yx)  , . . . ,  <y„) ,  (T))  .  (3) 

The  above  closure  fails  due  to  the  size  of  instantaneous  local  deviations  from  the  modelled 
unconditional  mean  values  in  the  flow  that  is  being  modelled.  The  CMC  model  seeks  to 
reduce  the  size  of  these  deviations  and  thereby  allow  a  similar  first  order  closure  (using 
different  averaging  operator)  to  succeed. 

The  CMC  model  seeks  to  employ  an  averaging  process  which  accounts  for  turbulence- 
induced  fluctuations.  Mixture  fraction  (f )  provides  a  useful  measure  of  the  degree  of  mixing 
between  fuel  and  oxidizer  masses,  which  is  invariant  under  chemical  reaction.  Averaging 
based  on  conditions  on  the  value  of  mixture  fraction  take  account  of  turbulence-induced 
mixing  fluctuations.  Mixture  fraction  is  a  scalar  which  is  derived  from  an  appropriate 
linear  combination  of  chemical  scalars, 

N 

t  =  E  (4) 

a—1 

such  that  the  net  chemical  production  of  that  scalar  is  identically  zero, 

n 

E  a«^cr  =  0  ,  (5) 

a=l 

and  is  normalised  so  that  (by  convention)  it  has  a  value  of  unity  in  fluid  that  originated 
wholly  from  the  fuel  stream,  and  zero  in  fluid  that  originated  in  the  oxidizer  stream.  It 
is  necessary  to  refer  to  fuel  and  oxidizer  ’origins’  since  while  fuel  and  oxidizer  may  be 
chemically  consumed  in  an  isolated  sample,  the  mixture  fraction  will  remain  unchanged. 
Thus  the  resultant  mixture  of  products  and  remnant  reactants  will  constitute  the  same 
value  of  mixture  fraction  as  the  original  unreacted  sample.  Only  the  introduction  of 
further  mass  (through  mixing)  to  change  the  atomic  population  of  the  sample  will  have 
any  bearing  on  mixture  fraction.  It  is  this  invariance  under  reaction  which  makes  mixture 
fraction  such  a  useful  coordinate  in  which  to  examine  nonpremixed  combustion. 
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The  CMC  method  involves  solving  for  the  mean  species  mass  fractions  which  result 
from  conditionally  averaging  an  ensemble  of  instantaneous  local  realisations  of  the  turbu¬ 
lent  flow  in  question.  The  conditional  mean  quantities  are  determined  by  discriminating 
between  contributing  samples  on  the  basis  of  their  associated  value  of  mixture  fraction. 
For  a  given  conditioning  value  (77)  of  mixture  fraction  (£),  only  those  samples  simultane¬ 
ously  meeting  the  conditional  criteria  (£(x,  t)  —  77)  are  allowed  to  contribute  to  the  mean. 
Thus  instead  of  traditional  uTxonditional  averaging,  where  averaging  over  an  ensemble 
of  points  in  time  and  space  would  only  produce  a  single  value,  the  conditional  averaging 
process  produces  a  profile  of  mean  values, 

Qa(x,t,rj)  =  (Ya{x,t)  1 77),  (6) 

as  a  function  of  the  conditioning  variable  abscissa  (77) . 

Conditional  averaging  of  the  instantaneous  local  species  mass  fraction  equations  (eq 
1)  yields, 

+  =  +  (/>>»» 1 7)  ,  P) 

where  Njj  is  the  conditional  mean  scalar  dissipation  rate  defined  by, 

ir>)  ■  (8) 

The  terms  in  equation  7,  from  right  to  left,  account  for  solution  unsteadiness,  mean  convec¬ 
tive  transport  in  physical  space,  turbulent  mixing  in  mixture  fraction  space,  and  chemical 
reaction.  Unlike  in  the  case  of  unconditionally  averaged  equations,  first  order  closure  of 
the  chemical  reaction  source  term  can  usually  be  achieved  by, 

(pwa(y1,...,yn,T)|77)«(p|77>d7a({y1i77),...,(yn|77),(T|77))  ,  (9) 

owing  to  the  smaller  degree  of  individual  sample  deviations  from  the  modelled  mean  quan¬ 
tities. 

In  practice  the  conditional  mean  scalar  dissipation  rate  must  be  determined  from  the 
conservation  equation  for  the  mixture  fraction  probability  density  function  (Pv), 

+  =  •  (io) 

to  ensure  conservation  of  mass  in  physical  space.  This  is  typically  done  by  using  an 
assumed  form  of  the  probability  density  function  (PDF)  such  as  a  beta  function  or  clipped 
Gaussian  distribution.  Equations  for  unconditional  mean  mixture  fraction  and  mixture 
fraction  variance  are  solved  in  physical  space  and  PDF  evolution  is  thus  determined  by 
the  evolution  of  these  two  controlling  parameters.  The  evolution  of  the  PDF  then  gives 
the  correct  form  of  the  conditional  mean  scalar  dissipation  rate, 

^  ■  (ll) 

It  can  be  seen  from  equation  8  that  Nv  is  a  non-negative  definite  distribution  in  mixture 
fraction  space.  The  boundary  conditions  upon  equation  11  are  well  known  [3,  9].  The 
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product  ( NVPV )  must  have  zero  value  and  zero  gradient  at  extreme  values  of  mixture 
fraction. 

Unconditional  mean  data  can  be  recovered  from  conditional  mean  profiles  through 
convolution  with  the  appropriate  mixture  fraction  PDF, 

(Ya(x,t))=  f  Pr,{x,t)(Ya{x,t)  |  T])dT)  .  (12) 

Jo 

2.1.1  Similarity  zones  in  conditional  statistics 

Past  applications  of  the  CMC  method  have  involved  the  use  of  assumptions  of  condi¬ 
tional  statistical  similarity  in  one  or  more  physical  directions  to  reduce  the  dimensionality 
of  the  problem  being  solved.  Steady  axisymmetric  turbulent  jet  flames  have  been  modelled 
by  assuming  that  conditional  statistics  are  invariant  across  the  radius  of  the  jet  at  any 
given  axial  location.  This  radial  independence  means  that  a  single  set  of  conditional  mean 
profiles  (one  for  each  species)  suffices  to  describe  the  behaviour  at  all  radial  points  at  a 
common  axial  location. 

Further,  due  to  the  parabolic  nature  of  the  flow  characteristics,  it  was  possible  to  use 
simple  marching  schemes  to  advance  the  solutions  from  nozzle  exit  plane  to  end  of  the 
flame  in  each  case.  In  effect,  these  simplifications  allowed  the  whole  steady-state  reactive 
flow  field  to  be  computed  progressively  using  a  single  cross-sectional  set  of  conditional 
mean  species  profiles,  which  is  marched  along  the  jet  axis  subject  to  the  CMC  equations. 

More  generally,  however,  these  types  of  assumptions  on  the  dimensional  degeneracy  of 
the  conditional  statistics  may  not  be  applicable.  In  a  worst  case  scenario,  an  individual 
conditional  mean  profile  would  have  to  be  solved  for  each  species  at  every  grid  location 
in  physical  space.  Typically,  between  thirty  and  one  hundred  grid  points  are  required  in 
each  mixture  fraction  profile  to  provide  adequate  resolution  of  nonpremixed  flame  struc¬ 
ture.  Thus  in  this  worst  case  scenario,  the  CMC  model  may  require  up  to  two  orders  of 
magnitude  more  computer  memory  and  very  much  more  computation  time  than  for  the 
solution  of  a  passive  scalar. 

For  problems  of  practical  interest,  such  as  gas  turbine  combustors,  it  is  important 
to  take  advantge  of  spatial  similarities  in  conditional  statistics  as  much  as  possible  to 
minimise  the  number  of  independent  profiles  which  must  be  carried  by  the  model.  It  is 
useful  to  introduce  the  concept  of  conditional  statistical  similarity  zones.  A  similarity 
zone  is  defined  here  as  a  region  of  space  across  which  conditional  statistics  are  spatially 
uniform  to  good  approximation.  Within  a  similarity  zone,  a  single  conditional  mean  profile 
for  each  species  suffices  to  describe  behaviour  at  all  physical-space  grid  points. 

In  Figure  2,  different  notional  boundaries  for  similarity  zones  are  depicted  for  a  variety 
of  flow  cases.  On  the  left  of  the  figure,  the  case  for  jet  flames  (described  above)  is  depicted, 
with  a  radial  marching  similarity  zone.  In  the  centre  of  the  figure,  possible  static  similarity 
zones  are  depicted  for  an  axisymmetric  combustor  can.  The  right  hand  side  of  Fig.  2  relates 
to  the  similarity  zone  approach  used  in  the  Conditional  Source  Evaluation  (CSE)  method 
[10],  which  is  described  in  Section  2.1.3. 

The  a  priori  determination  of  static  similarity  zone  boundaries  in  complex  flow  is 
arbitrary  and  somewhat  problematic.  In  two  dimensional  calculations,  streamlines  may 
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Simple  jet  Swirl  combustor  Swirl  combustor 

(CMC)  (CMC)  (CSE) 

Marching  radial  similarity  zone  Static  irregular  similarity  zones  Floating  similarity  zones 

Figure  2:  Schematic  diagram  indicating  the  differences  in  conditional  statistical  similarity 
zone  structure  in  jet  flames  and  swirl  combustors  for  CMC  and  CSE  methods. 


be  used  as  a  guide  in  the  apportioning  of  the  domain  into  similarity  zones.  As  with  the 
simple  jet  flow,  it  is  appropriate  that  where  streamlines  are  close  and  parallel,  a  similarity 
zone  can  be  delineated  across  these  streamlines.  Thus  in  a  locally-free  shear  layer,  zone 
boundaries  will  be  perpendicular  to  groups  of  parallel  streamlines.  In  order  to  capture 
conditional  mean  species  development,  zones  should  be  established  at  relatively  small 
intervals  along  streamlines.  Near  walls  and  other  flow  boundaries,  additional  similarity 
zones  may  be  required  to  capture  gradients  in  reactive  species.  The  establishment  of  a 
static  zonal  structure  will  probably  require  detailed  knowledge  of  the  flow  field  beforehand, 
perhaps  using  a  cold  flow  solution  as  a  guide. 

Given  the  definition  of  a  similarity  zone,  it  is  possible  to  average  a  divergence  form 
of  the  conditional  mean  species  equation  (eqn  7)  over  the  volume  (Vh)  of  such  a  zone  so 
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that, 


({ pui  1 77 )PvQq)  -  Qa 


d_ 

dxi 


{{pm  |  Tj)  Pv) 


NvPn 


d2Qc 
dr 12 


+  {pWa  |  77)  Pr, 


(13) 

where  {. .  ,}n  denotes  the  volume  averaging  operator.  It  is  possible  to  further  simplify  the 
above  equation  to  eliminate  one  of  the  left  hand  side  convective  terms.  However,  in  the 
light  of  the  potentially  irregular  boundaries  of  similarity  zones,  the  above  divergence  form 
is  easier  to  accurately  implement. 


Application  of  Gauss’s  Theorem  and  the  independence  of  conditional  mean  species 
data  upon  physical  location  give, 


(Fq  1 1)>-0„  (FPb)  =  {«■„}„ 


gg. 

dif 


+  (pWa  |  7?)  , 


(14) 


where  {Fq  \  rj)  and  ( Fp  \  rj)  denote  conditional  mean  species  mass  flux  and  conditional 
mean  mass  flux,  across  the  zone-bounding  surface.  These  fluxes  are  defined  as, 

(Fq  |  n)  =  jpFrr  js  (w  I  -)>'  KQ'Jsi  (is) 

(Fp  I  n)  =  T7TT- v-  f  1 ^ds‘  •  <16> 

where  Vh  is  the  volume  of  the  similarity  zone.  The  conditional  mean  quantities  at  the 
bounding  interface  have  intermediate  values  between  those  of  the  adjacent  zone  and  the 
current  zone.  Intermediate  valued  quantities  are  indicated  in  the  above  definitions  by 
superscript  prime  symbols.  The  companion  mixture  fraction  PDF  equation  can  be  treated 
similarly  to  the  CMC  equation  to  yield, 

=  •  (17> 


The  similarity  zone  formulation  embodied  by  equations  14-17  provides  a  means  for 
employing  the  CMC  method  in  cases  of  practical  interest.  Care  must  be  taken  in  imple¬ 
mentation  of  these  equations  to  maintain  their  adjoint  relationship.  Failure  to  do  so  would 
result  in  a  failure  to  conserve  mass  within  the  model. 


2.1.2  Determining  conditional  mean  auxiliary  variables 

Conditional  mean  quantities  such  as  velocity  {{m  \  v))  and  scalar  dissipation  rate  {Nv) 
must  be  determined  from  unconditional  mean  data  in  order  to  solve  equations  7-17.  Within 
a  conditional  statistical  similarity  zone  (see  Section  2.1.1)  these  variables  are  related  to 
the  unconditional  mean  variables  through  convolution  with  the  mixture  fraction  PDF, 

(u(z))  =  f  {u\t})  Pv(x)dri  +  eu(x)  (18) 

{N{x))=  f  Nr,Pj,{x)dTj+  ejv(z)  ,  (19) 

Jo 
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where  the  PDF  is  itself  usually  an  assumed-form  distribution,  dependent  upon  mixture 
fraction  mean  (£)  (x)  and  variance  (£/2  (£.))•  It  is  appropriate  to  include  the  small  error 
terms  eu(x)  and  ejv(x)  in  the  above  expressions  to  account  for  imperfect  conditional  sta¬ 
tistical  similarity  within  a  designated  zone.  In  order  to  determine  the  conditional  mean 
data,  a  solution  must  be  found  for  the  above  Fredholm,  integral  equations  of  the  first  kind 
[11]. 

Where  there  are  many  grid  points  in  a  similarity  zone  with  different  unconditional 
mean  values,  it  is  possible  to  employ  zeroth-order  linear  regularisation  [11]  to  estimate  the 
conditional  mean  data  for  the  zone  which  correspond  to  these  unconditional  mean  values. 
The  basis  for  linear  regularisation  is  perhaps  better  understood  when  a  discretised  form 
of  the  above  equations  is  written  as  a  linear  problem, 

A  •  u  +  e  =  b  ,  (20) 

where  the  matrix  A  represents  the  discretised  form  of  the  PDF  convolution  at  a  variety  of 
points  in  physical  space,  b  represents  the  unconditional  mean  values  at  points  in  physical 
space,  e  denotes  the  small  errors  at  points  in  physical  space  due  to  imperfect  conditional 
similarity,  and  u  represents  the  conditional  mean  values  in  mixture  fraction  space.  The 
matrix  A  is  usually  singular  since  the  convolution  is  essentially  a  smoothing  process  and, 
further,  their  is  no  requirement  for  A  to  be  square  in  any  case.  Thus  the  problem  is  likely 
to  be  ill-posed. 

Linear  regularisation  involves  forming  a  functional  of  u,  which  is  a  global  measure  of 
the  size  of  the  unknown  zonal  similarity  errors, 

/u  =  I  A  •  u  —  b  |2  (21) 

and  minimising  the  sum  of  this  functional  and  an  a  priori  functional, 

=  i  B  •  u  |2  (22) 

which  represents  a  global  measure  of  the  tendency  for  u  to  display  a  certain  profile  shape 
in  mixture  fraction  space.  The  composition  of  the  matrix  B  is  selected  according  to  the 
anticipated  form  of  the  conditional  mean  data,  be  it  essentially  constant,  linear,  parabolic, 
or  any  hybrid  of  these  forms  in  mixture  fraction  space  [11].  The  addition  of  the  a  priori 
information  places  further  constraints  on  the  form  of  u  so  that  the  linear  system  is  no 
longer  singular  and  a  unique  solution  is  then  possible  for  a  given  combination  of  the  two 
functionals.  The  overall  minimisation  of  fu  +  jgu,  subject  to  a  threshold  value  of  /u,  can 
be  written  as  a  normal  linear  equation  of  the  form, 

(At  •  A  -(-  7Bt  •  B)  •  u  =  At  •  b  (23) 

where  7  is  the  blending  factor  with  0  <  7  <  00  through  which  the  minimisation  is  achieved. 

It  can  be  seen  that  very  large  7  values  will  cause  minimisation  to  be  based  solely  on  u 
conforming  to  the  a  priori  assumed  shape,  rather  than  matching  the  unconditional  mean 
data.  Conversely,  very  small  7  values  cause  minimisation  to  be  based  solely  on  matching 
the  unconditional  mean  data,  which  as  stated  above,  will  not  have  a  unique  solution  if 
the  number  of  points  in  mixture  fraction  space  exceeds  the  number  of  physical  grid  points 
present  in  the  zone.  Where  a  unique  solution  does  exist  in  cases  of  7  — >■  0,  it  often  exhibits 
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wild  oscillations  that,  while  satisfying  the  unconditional  mean  data  set,  are  not  physically 
plausible. 

Press  et  al  [11]  suggest  that  the  blending  factor  should  be  selected  so  that  the  lower 
threshold  values  of  functional  fu  should  not  be  less  than  an  estimate  of  the  overall  error 
associated  with  the  eu(x)  terms.  To  have  an  unreasonably  lower  threshold  value  invites 
unphysical  profile  behaviour  as  a  result  of  too  many  degrees  of  freedom  in  the  conditional 
mean  profile.  Estimates  of  the  overall  error  which  should  be  associated  with  eu  and  e/v 
can  be  made  based  on  the  unconditional  mean  values  for  each  variable  and  the  coarseness 
of  the  similarity  zones  employed. 

There  are  two  alternate  means  by  which  the  CMC  and  PDF  equations  (14,  17)  can  be 
solved  using  the  auxiliary  variables  ( U{  \  if)  and  Nv.  The  first  involves  determining  (u;  |  rj) 
from  linear  regularisation  of  the  unconditional  mean  velocity  fields  ((u,-))  within  each 
similarity  zone.  The  conditional  mean  velocity  is  then  used  in  conjunction  with  equation 
17  to  find  the  conditional  mean  scalar  dissipation  rate  {Nj,}^).  These  two  quantities  are 
then  used  in  the  solution  of  equation  14.  This  procedure  is  preferrable  since  it  retains 
the  strong  link  between  the  two  equations,  required  to  conserve  mass.  This  procedure  can 
cause  some  difficulty  in  ensuring  boundary  conditions  and  realisability  conditions  are  met 
for  Nv.  Smith  [3]  describes  some  of  the  difficulty  which  can  arise  from  solving  for  Nv  from 
the  PDF  equation  directly.  It  was  found  that  in  many  cases,  inaccurate  solutions  for  the 
(£)  and  (£'2)  fields,  which  determine  Pv  and  its  spatial  variation,  can  generate  unrealisable 
Nv  profiles  which  can  quickly  destabilise  the  CMC  equations.  These  inaccurate  solutions 
are  likely  in  the  early  iterations  of  an  iterative  steady  state  solver,  and  some  means  must 
be  available  to  otherwise  approximate  Nn  under  these  circumstances. 

The  second  means  of  determining  Nv  is  far  less  prone  to  failing  the  realisability  con¬ 
straints.  This  procedure  involves  determining  Nv  by  linear  regularisation  from  the  un¬ 
conditional  mean  scalar  dissipation  rate  ( N )  which  is  generated  in  the  solution  of  the 
(f'2)  equation.  This  second  technique  can  provide  provisional  Nv  profiles  for  use  in  the 
CMC  equations  where  the  first  technique  is  yielding  non-physical  profiles  due  to  poorly 
converged  mixture  fraction  fields. 

2.1.3  CSE  method 

The  Conditional  Source  Evaluation  (CSE)  method  was  proposed  by  Bushe  and  Steiner 
[10]  as  a  use  of  linear  regularisation  to  predict  chemical  reaction  rates  in  turbulent  non- 
premixed  combustion.  The  CSE  method  seeks  to  remedy  the  classical  closure  problem  (see 
equation  3)  directly  without  having  to  solve  the  CMC  equations.  Instead,  chemical  species 
transport  and  mixing  is  solved  in  terms  of  unconditional  mean  variables  (as  in  equation  2). 
At  each  step  where  unconditional  mean  chemical  reaction  rates  are  required,  conditional 
mean  species  profiles  are  determined  from  the  unconditional  mean  species  data  using  lin¬ 
ear  regularisation,  and  this  conditional  mean  information  is  used  to  evaluate  conditional 
mean  reaction  rates.  These  reaction  rates  are  convolved  with  mixture  fraction  PDFs  to 
return  unconditional  mean  reaction  rates. 

The  success  or  failure  of  the  method  depends  to  a  very  large  extent  on  the  accuracy  of 
the  linear  regularisation  which  determines  the  conditional  mean  species  profiles.  As  with 
the  similarity  zone  formalism  in  CMC  modelling,  CSE  also  requires  a  similarity  zone  in 
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order  to  allow  linear  regularisation  to  take  place.  Static  similarity  zones  (see  Fig.  2)  can 
be  employed  as  in  the  CMC  model  with  minimal  computational  effort.  Alternatively  the 
CSE  method  can  employ  a  floating  similarity  zone  for  each  physical  grid  point.  Each  of 
these  similarity  zones  then  includes  the  local  grid  points  surrounding  the  point  for  which 
reaction  rates  are  being  determined. 

The  major  problem  with  the  CSE  method  is  that  the  regularisation  process  can  produce 
conditional  mean  species  profiles  which  violate  atomic  mass  balances,  resulting  in  impos¬ 
sible  chemical  compositions.  These  compositions  can  destabilise  the  system  of  chemical 
reactions  causing  the  computation  to  diverge. 

A  simple  test  can  be  derived  to  determine  the  accuracy  of  the  CSE  method.  By 
taking  known  conditional  mean  chemical  profiles  (say  resulting  from  a  chemical  equilibrium 
calculation),  and  convolving  them  with  known  PDFs,  a  set  of  unconditional  mean  data  can 
be  derived.  The  test  for  the  CSE  method  is  to  return  the  conditional  mean  data  from  the 
unconditional  mean  data.  Such  a  test  was  performed  using  species  profiles  corresponding 
to  equilibrium  conditions  for  C2H4  and  air  at  a  variety  of  equivalence  ratios. 


R(mj  Rfa] 

Figure  3:  Plots  of  mixture  fraction  unconditional  mean  and  variance  as  functions  of  radial 
position  at  various  axial  locations  in  a  30  m/s  C2H4  jet  flame.  Symbols  denote  axial 
location  relative  to  the  stoichiometric  flamelength  as  :  +  -  x/L  =  1/3,  X  -  x/L  —  2/3  , 
and  *  -  x/L  =  1 

The  mixture  fraction  PDF  values  were  determined  to  be  beta  functions  depending  upon 
the  mean  and  variance  of  mixture  fraction  at  a  series  of  physical  points  taken  from  radial 
slices  through  an  axisymmetric  fast  chemistry  jet  flame  calculation.  Figure  3  contains 
plots  of  the  unconditional  mean  and  variance  of  mixture  fraction  as  functions  of  the  radial 
location  in  each  of  these  slices.  It  can  be  seen  that  the  centreline  mean  mixture  fraction 
decays  with  increasing  axial  distance  from  the  nozzle,  while  the  half-width  of  the  jet 
spreads,  as  jet  fluid  is  mixed  with  surrounding  fluid.  The  mixture  fraction  variance  peaks 
upstream  near  the  interface  between  the  pure  fuel  and  pure  oxidizer  streams  in  the  jet 
flow.  The  variance  rapidly  decays  as  a  result  of  diffusive  dissipation  as  the  fluid  moves 
downstream. 

Using  chemical  equilibrium  conditional  mean  profiles,  and  the  mixture  fraction  PDFs 
resulting  from  the  above  mixture  fraction  distributions,  unconditional  mean  species  data 
was  generated  as  input  to  the  CSE  integral  equation  solver.  The  a  priori  information 
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Figure  f:  Temperature  and  OH  radical  mass  fraction  as  functions  of  mixture  fraction  at 
various  axial  locations  in  a  C2H4  jet  flame.  From  the  left  to  right ,  the  plots  correspond  to 
x/L  =  1/3, 2/3, 1.  In  each  figure,  the  symbols  denote,  x  a  priori  conditional  mean  data, 
+  CSE  model  estimate  of  a.  priori  data,  and  *  unconditional  mean  species  data  determined 
from  a  priori  data  and  used  as  input  to  CSE. 


matrix  B  was  carefully  constructed  for  each  species  to  select  for  the  correct  profile  shape 
in  mixture  fraction  space.  For  regions  where  the  equilibrium  conditional  mean  data  was 
essentially  constant,  that  portion  of  the  B  matrix  was  made  to  represent  a  first  derivative 
operation.  Regions  of  essentially  linear  behaviour  were  represented  by  a  second  deriva¬ 
tive  operator,  and  regions  of  parabolic  behaviour  were  represented  by  a  third  derivative 
operator.  Minimisation  subject  to  B  thus  sought  to  minimise  slope  in  regions  thought  to 
be  constant,  curvature  in  regions  thought  to  be  linear,  and  change  in  curvature  in  regions 
thought  to  be  parabolic. 

Figure  4  shows  the  results  of  the  CSE  integral  equations  compared  with  the  a  pri¬ 
ori  conditional  mean  profiles  which  they  should  match  and  the  input  unconditional  mean 
profiles  generated  by  the  a  priori  conditional  mean  information.  The  two  reactive  scalars 
plotted  are  indicative  of  the  type  of  agreement  found  for  all  of  the  chemical  species  em¬ 
ployed.  In  each  plot,  the  range  in  mixture  fraction  plotted  corresponds  with  the  range 
in  mean  mixture  fraction  found  in  the  radial  slice  in  question.  As  no  unconditional  data 
exists  for  mixture  fractions  outside  this  range  it  is  not  possible  to  accurately  determine 
conditional  mean  profiles  in  these  regions  using  linear  regularisation. 

It  is  apparent  that  the  agreement  between  the  imposed  and  modelled  conditional  mean 
temperature  profiles  is  good  at  all  axial  locations  but  improves  somewhat  at  greater  dis¬ 
tances  from  the  nozzle.  The  agreement  is  poorer  for  the  radical  species  OH,  with  the 
modelled  profiles  containing  unphysical  negative  mass  fractions  at  some  mixture  fractions. 
The  agreement  between  the  OH  profiles  does,  however,  improve  at  the  downstream  lo¬ 
cations.  The  reason  for  the  gradual  improvement  in  agreement  for  all  species  with  axial 
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distance  from  the  nozzle,  is  that  the  mixture  fraction  variance  decays  so  that  uncondi¬ 
tional  mean  values  are  more  representative  of  conditional  mean  values.  Integral  equation 
solutions  are  most  easily  found  when  the  unconditional  mean  data  points  are  separated  in 
mixture  fraction  space  by  many  standard  deviations  of  mixture  fraction.  Lower  variances 
give  rise  to  a  smaller  degree  of  overlap  between  the  mixture  fraction  PDFs  corresponding 
to  different  physical  locations. 

The  unphysical  appearance  of  negative  mass  fractions  requires  special  handling  before 
using  the  conditional  mean  species  profiles  to  determine  conditional  mean  reaction  rates. 
Simple  clipping  of  negative  mass  fractions  to  zero  may  be  sufficient  to  avoid  numerical 
instability.  The  results  shown  above  indicate  that  the  CSE  method  may  be  a  useful 
alternative  to  the  CMC  model  in  the  prediction  of  turbulent  nonpremixed  combustion. 


2.2  Heat  Transfer 


To  date,  CMC  modelling  has  been  limited  to  applications  of  free-shear  flows,  such  as 
turbulent  jet  flames.  The  absence  of  walls  in  these  applications  obviates  the  need  for  a 
means  of  treating  conductive  heat  transfer  to  solid  surfaces.  However,  practical  application 
of  CMC  modelling  demand  a  capability  for  dealing  with  heat  transfer  to  adjacent  walls. 

In  the  immediate  vicinity  of  a  wall,  some  of  the  underpinning  assumptions  of  the  CMC 
method  break  down.  The  existence  of  a  laminar  sub-layer  preclude  the  high  Reynolds 
number  turbulent  transport  approximations  which  permeate  the  model.  Thus  conditional 
similarity  zones  cannot  be  extended  to  walls.  There  must  be  an  interface  between  the  wall 
and  the  adjacent  similarity  zone  boundary.  Heat  diffuses* across  this  interface  to  the  wall 
according  to, 

<*>-!§?•  <24) 


where  (sc)  denotes  the  unconditional  mean  energy  loss  from  the  fluid  due  to  conduction, 
and  (gf)  denotes  the  unconditional  mean  heat  flux  in  the  i-th  direction  due  to  conduction. 
The  heat  flux  is  given  by, 

«>  =  >  <25> 

where  A  is  the  thermal  conductivity,  a  function  of  molecular  activity,  and  (T)  is  the 
unconditional  mean  temperature. 

Note  that  in  the  purely  turbulent  portion  of  the  flow,  heat  conduction  by  molecular 
transfer  is  effectively  incorporated  into  a  turbulent  thermal  diffusivity  term,  which  has  a 
conditional  mean  counterpart  in  the  CMC  equations  for  standardised  enthalpy.  At  the 
edge  of  turbulent  region  special  boundary  terms  must  be  employed  to  account  for  this 
laminar  heat  transfer  effect.  It  is  proposed  that  a  conditional  mean  conductive  loss  term 
{<*  |  rj)}n  be  employed  in  wail-boundary  zones. 

In  cases  with  specified  wall  temperatures,  the  unconditional  mean  heat  loss  across 
the  wall  interface  can  be  determined  using  unconditional  mean  temperatures  at  the  zone 
boundary.  These  temperatures  can  be  determined  using  the  convolution, 


(T(£)>=  [l{{T\V)}aP„(lL)dr,  . 
Jo 


13 


DSTO-TN-0256 


The  computed  unconditional  mean  sources  along  the  interface  can  then  be  employed  to 
determine  conditional  mean  sources  using  linear  regularisation  (see  Section  2.1.2)  of  the 
following, 

(sc{x))=  f  {(sc\T))}uPv{x)dT)  .  (27) 

Jo 

In  the  case  of  wall  boundaries  where  the  heat  flux  is  specified,  the  wall  temperatures  are 
inferred  from  the  unconditional  mean  temperatures  on  the  adjacent  bound  of  the  adjacent 
similarity  zone.  Regularised  solutions  for  conditional  mean  heat  flux  are  applied  to  the 
CMC  equations  as  in  the  specified  temperature  boundary  condition. 

It  is  worth  noting  that  it  is  also  possible  to  use  this  procedure  to  account  for  unusual 
effects  such  as  chemical  species  migrating  to  and  from  the  wall  as  a  result  of  molecular 
diffusion.  Such  an  approach  will  be  useful  if  surface  chemistry  (like  platinum-film  catalysis) 
is  to  be  treated.  Provided  the  overall  mass  balance  between  the  wall  and  turbulent  fluid  is 
maintained  this  process  will  not  influence  the  conserved  scalar  field  upon  which  conditional 
averaging  is  predicated. 

2.2.1  Radiation  heat  transfer 

Radiation  heat  transfer  has  been  modelled  in  the  past  using  the  assumption  of  an 
optically  thin  emitting  media.  Under  this  assumption,  radiation  emitted  by  hot  gases  is 
transmitted  from  the  domain  without  any  reabsorption  or  scattering  whilst  in  transit. 
This  assumption  reduces  radiant  heat  loss  to  being  solely  a  function  of  local  conditions 
rather  than  a  non-local  function  of  conditions  throughout  the  domain. 

The  local  nature  of  the  optically  thin  radiation  treatment  allows  a  conditional  mean 
heat  loss  term  to  be  applied  in  the  CMC  equations  without  any  need  for  describing  geo¬ 
metric  disposition  of  the  domain.  Thus 

(sr  |  77)  «  47rcr  kp  ((T  ]  p) ,  (YH20  i  v)  i  (Yc02  I  v)  ?  (Y soot  I  v)  >  •  •  •)  (T  I  l)4  »  (2^) 

applies  where  <7  denotes  the  Stefan-Boltzmann  constant,  and  kp  denotes  the  Planck  mean 
absorption  coefficient  for  the  mixture  at  the  given  mixture  fraction. 

The  optically  thin  assumption  has  been  found  to  be  a  good  approximation  in  many 
non-sooting  combustion  applications  where  physical  dimensions  are  of  a  laboratory  or 
aero-engine  scale.  The  possibility  exists,  however,  that  the  assumption  will  break  down  in 
the  presence  of  soot  and  that  parts  of  the  domain  may  absorb  energy  transmitted  from 
other  parts  of  the  domain.  Under  these  circumstances,  the  geometric  disposition  of  the 
participating  media  and  the  domain  boundaries  must  be  included  in  the  calculation  of 
radiant  transfer. 

A  Finite  Volume  Radiation  (FVR)  model  [12,  13]  can  be  employed  under  these  cir¬ 
cumstances  to  predict  radiation  heat  transfer.  Flux  intensities  must  be  calculated  for  each 
incremental  solid  angle  of  direction  for  the  boundaries  of  each  grid  cell.  The  balance  of 
transmitted  intensity  into  and  out  of  each  direction  in  each  cell  is  made  up  by  the  net 
emission  or  absorption  within  that  cell.  The  result  is  a  large  system  of  unconditional  mean 
intensity  equations  which  must  be  inverted  at  each  step.  Though  tedious,  this  process  is 
relatively  straight-forward. 
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AXIAL  LOCATION  [m] 


Figure  5:  Comparison  of  unconditional  mean  centreline  temperature  variations  with  axial 
distance  from  the  nozzle  in  a  sooting  C2H4  jet  flame,  for  different  radaition  heat  loss 
models.  Symbols  denote  models  according  to  :  +  no  model,  x  optically-thin  model, 

and  □  FVR  model. 

Determination  of  the  unconditional  mean  fourth-power  of  temperature  is  required  in 
each  grid  cell  in  order  to  determine  unconditional  mean  absorption  and  emission.  This 
can  be  determined  from  the  convolution  of  the  conditional  mean  temperature,  raised  to 
the  fourth-power,  with  the  mixture  fraction  PDF  in  each  grid  cell.  Conversely,  the  energy 
source/sink  effect  upon  the  CMC  system  of  equations  is  that  the  unconditional  mean 
emission  or  absorption  in  each  grid  cell  must  be  regularised  to  generate  conditional  mean 
emission/absorption  terms.  The  regularisation  follows  the  same  process  as  that  described 
in  the  preceding  sections. 

While  feasible,  use  of  radiation  models  of  the  complexity  of  the  FVR  models,  adds  a 
large  computational  burden  to  the  already  expensive  task  of  turbulent  combustion  mod¬ 
elling.  It  is  worthwhile  in  most  cases  to  determine  ahead  of  time  if  use  of  these  models  is 
warranted  over  a  vastly  simpler  optically-thin  radiation  model.  This  can  be  done  using  a 
simple  fast  chemistry  or  laminar  flamelet  model  for  turbulent  combustion  (see  [7]).  These 
models  have  poor  predictive  capability  in  most  problems  of  interest,  but  are  computation¬ 
ally  cheap  and  will  give  good  order-of-magnitude  estimates  for  the  purposes  of  determining 
if  an  FVR  model  is  required  for  the  treatment  of  radiation. 

This  type  of  analysis  has  been  employed  in  the  following  example  for  a  turbulent 
axisymmetric  C2H4  sooting  jet  flame  surrounded  by  adiabatic  walls  at  a  radial  distance 
of  80  nozzle  diameters.  Laminar  flamelet  calculations  were  made  where  radiant  heat  loss 
to  the  walls  was  computed  using  optically-thin  and  FVR  treatments.  Figure  5  provides 
a  comparison  of  the  centreline  unconditional  mean  temperature  ((T))  distributions  in  the 
flame  which  result  from  the  application  of  the  FVR  model,  the  optically-thin  model,  and 
no  model  at  all  (no  radiation  heat  loss).  It  is  apparent  that  the  difference  between  the 
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results  predicted  by  the  two  radiation  models  is  very  slight  in  this  case  compared  to  the 
results  for  an  adiabatic  flame.  The  short  path  lengths  for  radiation  emission  through  the 
flame  in  this  computation,  and  the  uniform  bounding  geometry  make  the  optically-thin 
model  a  good  choice  for  such  a  case. 


2.3  Chemical  Reactions 

The  presence  of  chemical  reaction  rates  in  the  CMC  equations  (Eq  13),  with  timescales 
which  vary  by  many  orders  of  magnitude,  cause  these  equations  to  be  numerically  stiff. 
Stiff  equation  sets  are  those  where  there  is  a  wide  disparity  between  the  largest  and 
smallest  timescales  associated  with  source  terms.  The  scales  of  interest  in  most  systems  are 
those  which  change  relatively  slowly  and  iterative  solution  steps  have  a  size  in  accordance 
with  these  timescales.  Stiff  systems  of  equations  require  the  iterative  step  sizes  to  be 
commensurate  with  the  shortest  timescales,  and  yet  the  duration  of  the  calculation  must 
still  run  for  a  period  in  step  with  the  longest  scales. 

Special  solution  algorithms  must  be  employed  to  solve  numerically  stiff  systems  of 
equations  which  require  a  great  deal  more  computational  effort  than  the  simple  solvers 
which  are  employed  for  pure  fluid  dynamic  problems.  In  practice,  direct  integration  of 
these  source  terms  typically  requires  95  —  99%  of  the  overall  CPU  time  devoted  to  an 
entire  computation  including  Favre-averaged  turbulence  modelling,  and  the  other  terms 
in  the  CMC  equations.  Stiffness  can  be  largely  avoided  in  cases  which  employ  reduced 
chemical  reaction  mechanisms.  Those  reduced  schemes  which  eliminate  the  computation 
of  key  radical  species  such  as  H,  OH  and  OH  also  eliminate  the  fastest  reactions  which 
are  the  cause  of  worst  of  the  stiffness.  However,  in  eliminating  these  species  from  direct 
computation,  their  important  influence  in  determining  many  combustion  characteristics  is 
lost  from  the  model. 

Another  means  of  decreasing  the  amount  of  CPU  time  spent  in  integrating  chemical 
source  terms  involves  in  situ  adaptive  tabulation  (ISAT).  The  results  of  integration  of 
chemical  terms  are  stored  in  a  dynamic  tree  structure  as  the  computation  proceeds.  Each 
integration  is  stored  according  to  its  start-point  chemical  composition  and  incremental 
time  interval  of  integration.  Computations  are  only  made  if  a  search  of  the  tree  structure 
reveals  that  a  similar  computation  has  not  been  already  made.  If  a  similar  computation 
has  been  made,  then  the  end-point  chemical  compositions  of  that  computation  are  re-used. 
The  definition  of  similar  can  be  taken  to  mean  that  the  search  composition  is  within  a  short 
range  of  a  stored  start-point  composition;  short  enough  so  that  the  results  of  integration 
of  the  search  composition  would  result  in  end-point  composition  within  a  small  distance 
from  the  stored  end-point  conditions. 

A  tree  structure  is  used  for  data  storage  instead  of  an  indexed  table  due  to  the  excessive 
dimensionality  that  the  latter  would  entail.  An  indexed  table  requires  a  separate  array 
dimension  for  each  chemical  species  as  well  as  enthalpy  and  the  interval  of  integration. 
Owing  to  the  nature  of  the  chemical  reactions  and  the  proliferation  of  unrealisable  chemical 
states,  the  hyperspace  described  by  such  a  table  would  be  mostly  empty  and  thus  storage 
efficiency  would  be  low. 

The  tree  structure  instead  uses  characteristic  chemical  species,  enthalpy  and  time  inter¬ 
val  values  to  establish  a  test  function  as  the  sum  of  normalised  entries  in  the  composition 
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SEARCH  AMONG  NODES 
AND  LEAVES  FOR  CLOSEST 


•  NODE 


COMPUTE  NEW  SOLUTION, 

IS  ERROR  GREATER  THAN  TOLERANCE  ? 


Figure  6:  Graphic  representation  of  the  leaf  and  node  creation  algorithm  employed  in  the 
in  situ  adaptive  tabulation  scheme. 


/  enthalpy  /  time  vector.  The  searching  process  proceeds  by  computing  the  value  of  this 
test  function  for  the  search  composition.  This  value  is  then  compared  to  values  at  tree 
nodes  to  the  ’left’  and  ’right’  of  the  current  search  node.  A  node  is  a  storage  location 
where  values  of  the  test  function  at  nodes  or  leaves  to  the  ’left’  and  ’right’  are  stored. 
The  next  search  node  becomes  that  to  which  the  search  value  is  closest.  The  difference 
between  the  new  ’left’  and  ’right’  nodes  to  the  new  -search  node  is  then  smaller  than  the 
preceding  difference,  and  so  the  search  closes  in  on  the  closest  stored  composition  to  the 
search  composition.  Eventually  the  search  reaches  a  stage  where  ’leaves’  (accessible  stored 
data)  branch  from  nodes  and  one  of  these  leaves  is  selected  as  the  closest  data  to  the 
search  data. 

The  treatment  of  leaf  data  upon  finding  a  matching  leaf  is  illustrated  by  the  diagram 
in  Figure  6.  Once  a  leaf  has  been  found,  the  size  of  the  discrepancy  between  the  search 
composition  and  stored  composition  is  determined.  If  this  discrepancy  is  within  a  small 
range,  then  the  stored  end-point  conditions  are  used  as  the  integrated  solution.  If  the 
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discrepancy  is  too  large,  then  integration  of  the  search  composition  is  undertaken.  If  the 
end-point  composition  from  this  integration  is  dissimilar  from  the  end-point  composition 
of  the  leaf  data,  then  a  new  leaf  is  established  to  save  the  newly-integrated  result.  In  so 
doing,  a  new  node  is  established  between  the  old  node  and  old  leaf,  and  connectivity  is 
reorganised  so  that  the  old  node  links  to  the  new  node  instead  of  the  old  leaf,  and  the  new 
node  links  to  the  old  leaf  and  the  new  leaf.  If  the  found-leaf’s  end-point  composition  is 
similar  to  the  computed  end-point  composition,  then  no  new  leaf  is  generated  and  instead 
the  leaf’s  start-point  range  is  expanded  to  include  the  start-point  search  composition.  In 
this  instance  the  existing  tree  connectivity  is  retained. 

The  size  of  the  end-point  composition  error  tolerance  is  a  determining  factor  in  the 
number  of  leaves  which  will  be  generated  in  a  given  computation.  If  this  error  tolerance 
is  large,  then  less  leaves  will  be  generated  since  it  will  be  more  probable  that  search  com¬ 
positions  will  fall  within  existing  leaf  ranges.  Thus  the  overall  computation  will  proceed 
faster  since  less  integration  is  required,  but  at  the  expense  of  somewhat  lower  accuracy. 
Naturally,  the  reverse  is  true  in  that  narrower  error  tolerances  will  tend  to  admit  smaller 
errors,  but  cause  more  leaves  to  be  generated  and  thus  more  integration  steps  to  be  taken. 

The  amount  of  computer  memory  available  limits  the  size  to  which  the  tree  can  be 
allowed  to  grow.  By  recording  the  number  of  times  individual  leaves  have  been  accessed 
and  individual  nodes  have  been  traversed,  it  is  possible  to  determine  which  parts  of  the 
tree  are  little  used.  These  branches  can  be  pruned  as  memory  limits  are  reached  to  make 
way  for  new  data  produced  by  current  calculations. 

Trees  generated  from  prior  computations  with  similar  chemistry  can  be  used  in  new 
computations  instead  of  growing  trees  from  scratch.  The  amount  of  time  saved  in  using 
an  existing  tree  can  be  substantial.  If,  for  example,  a  computation  is  repeated  with  the 
same  chemistry  on  a  slightly  different  grid,  then  one  might  reasonably  expect  that  very 
few  new  chemical  integrations  will  be  required.  In  this  case,  computation  time  may  be 
reduced  by  orders  of  magnitude. 

For  nonpremixed  combustion  it  has  been  found  that  it  is  prudent  to  employ  a  number 
of  trees,  each  storing  data  for  chemical  reactions  occurring  at  widely  different  stoichiome¬ 
tries.  This  is  done  because  of  the  very  different  nature  of  fuel-rich  reactions  and  fuel-lean 
reactions  for  many  fuels.  To  store  the  entire  range  of  stoichiometries  in  a  single  tree 
can  result  in  unnecessarily  long  search  times.  Trees  can  be  indexed  according  to  mixture 
fraction  ranges,  and  tree  searching  can  proceed  within  those  broad  ranges. 

Past  CMC  computations  of  turbulent  jet  flames  using  ISAT  for  reaction  mechanisms 
with  32  species  and  111  reactions,  have  shown  up  to  a  ~  25%  saving  in  computation 
time.  It  should  be  pointed  out  however  that  these  types  of  computations  use  a  parabolic 
marching  similarity  zone  to  determine  a  steady  state  solution  (see  Section  2.1.1).  This 
marching  procedure  does  not  call  for  repeated  iterations  of  the  same  flow  field  and  is  thus 
not  likely  to  derive  much  benefit  from  ISAT.  It  is  expected  that  elliptic  calculations  of  the 
type  required  in  gas  turbine  combustor  computations  are  likely  to  be  accelerated  by  ISAT 
to  a  substantially  greater  degree. 
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3  Concluding  Remarks 

The  key  physical  processes  which  contribute  to  modelling  gas  turbine  combustors 
(GTCs)  have  been  outlined.  Taking  a  combustion  model  as  a  starting  point,  it  has  been 
described  how  the  effects  of  turbulent  mixing  and  heat  transfer  can  be  included. 

The  Conditional  Moment  Closure  (CMC)  model  for  turbulent  nonpremixed  combustion 
has  been  generalised  to  allow  for  its  application  in  a  GTC  configuration,  where  complex 
recirculating  flow,  high  turbulence  levels,  and  significant  heat  transfer  to  bounding  walls 
are  all  important.  Perhaps  the  key  aspect  of  this  generalisation  is  the  concept  of  conditional 
similarity  zones  which  allow  auxiliary  conditional  mean  quantities  to  be  computed  from 
conventionally  averaged  quantities,  using  linear  regularisation  of  the  associated  sets  of 
integral  equations. 

A  possible  simpler  alternative  to  the  CMC  model,  the  Conditional  Source  Evaluation 
(CSE)  method  has  been  described  and  evaluated.  The  CSE  method  makes  extensive  use 
of  linear  regularisation  to  solve  for  conditional  mean  chemical  species  profiles  directly.  The 
method  shows  some  promise  but  requires  further  evaluation. 

Finally,  it  was  recognised  that  the  numerical  stiffness  of  the  CMC  equations  can  be  sig¬ 
nificant  and  that  repeated  direct  itegration  of  such  stiff  equations  for  very  similar  chemical 
steps  is  wasteful  of  resources.  An  adaptive  tabulation  method  has  been  described  to  reduce 
the  cost  of  solving  chemical  reaction  steps  in  time,  by  storing  and  retrieving  computed 
steps  in  a  dynamic  tree-storage  structure. 

Other  aspects  of  GTC  modelling,  such  as  liquid  fuel  and  other  multi-phase  effects,  are 
also  under  current  development  [14,  15].  Those  aspects,  which  have  been  omitted  from 
discussion  here,  will  also  be  implemented  along  with  the  models  discussed  in  this  note,  in 
the  remainder  of  Task  DST  98/091. 
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